#To make "figure 4" and columns 6

rm(list=ls()) #shortcut to remove data
setwd("~/Dropbox/research/ATUS")

library(readstata13)
library(wesanderson)

all_atus <- read.dta13("all_atus.dta")
all_atus<-all_atus[!(all_atus$age<18 | all_atus$age>65 | all_atus$unclassified_paper>0),] #working population

#For the figures:
#Use trend1, trend2, and trend 3 to create the ablines in the figures
trend1<- lm(work_paper ~ year, data=subset(all_atus), year <=2008, weights=weight_adj) 
trend2<- lm(home_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
trend3<- lm(leisure_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 

#Take the mean of every year, and create a matrix for plot: #weighted mean.
m03 <- mean(all_atus$work_paper[all_atus$year==2003], w=all_atus$weight_adj[all_atus$year==2003])
m04 <- mean(all_atus$work_paper[all_atus$year==2004], w=all_atus$weight_adj[all_atus$year==2004])
m05 <- mean(all_atus$work_paper[all_atus$year==2005], w=all_atus$weight_adj[all_atus$year==2005])
m06 <- mean(all_atus$work_paper[all_atus$year==2006], w=all_atus$weight_adj[all_atus$year==2006])
m07 <- mean(all_atus$work_paper[all_atus$year==2007], w=all_atus$weight_adj[all_atus$year==2007])
m08 <- mean(all_atus$work_paper[all_atus$year==2008], w=all_atus$weight_adj[all_atus$year==2008])
m09 <- mean(all_atus$work_paper[all_atus$year==2009], w=all_atus$weight_adj[all_atus$year==2009])
m10 <- mean(all_atus$work_paper[all_atus$year==2010], w=all_atus$weight_adj[all_atus$year==2010])
m11 <- mean(all_atus$work_paper[all_atus$year==2011], w=all_atus$weight_adj[all_atus$year==2011])
m12 <- mean(all_atus$work_paper[all_atus$year==2012], w=all_atus$weight_adj[all_atus$year==2012])
m13 <- mean(all_atus$work_paper[all_atus$year==2013], w=all_atus$weight_adj[all_atus$year==2013])
m14 <- mean(all_atus$work_paper[all_atus$year==2014], w=all_atus$weight_adj[all_atus$year==2014])
mean <- cbind(m03,m04, m05, m06, m07, m08, m09, m10, m11,m12, m13, m14)

year <- seq(2003, 2014, 1)

#Format properly
mat<- rbind(year, mean) #bind mean to year
mat <- t(mat) #transpose matrix, so it's 2 rows of year and mean
colnames(mat)[colnames(mat)==""] <- "mean"

#diagram for work_paper
#plot(mat[,1], mat[,2], type = "o") #plot takes in the xs and ys and draws a line by type 0
                                   #xs are years, ys are trend.

mat.d <- as.data.frame(mat)

# market work trend plot ####
#pal <- wes_palette("Zissou", 100, type = "continuous")

ggplot(mat.d, aes(x=mat[,1], y=mat[,2])) +
       geom_point(shape=1) + 
       geom_smooth(method=lm, linetype = 2, se= FALSE) +
       geom_line(size = 1) +
#       scale_y_continuous(limits = c(20, 30)) +
       scale_x_continuous(breaks = 2003:2014) +
       xlab("Year") +
       ylab("Hours Per Week") +
       ggtitle("Market Work")

#Repeat for the other two variables to be diagramed:
#Take the mean of every year, and create a matrix for plot:
m03 <- mean(all_atus$home_paper[all_atus$year==2003], w=all_atus$weight_adj[all_atus$year==2003])
m04 <- mean(all_atus$home_paper[all_atus$year==2004], w=all_atus$weight_adj[all_atus$year==2004])
m05 <- mean(all_atus$home_paper[all_atus$year==2005], w=all_atus$weight_adj[all_atus$year==2005])
m06 <- mean(all_atus$home_paper[all_atus$year==2006], w=all_atus$weight_adj[all_atus$year==2006])
m07 <- mean(all_atus$home_paper[all_atus$year==2007], w=all_atus$weight_adj[all_atus$year==2007])
m08 <- mean(all_atus$home_paper[all_atus$year==2008], w=all_atus$weight_adj[all_atus$year==2008])
m09 <- mean(all_atus$home_paper[all_atus$year==2009], w=all_atus$weight_adj[all_atus$year==2009])
m10 <- mean(all_atus$home_paper[all_atus$year==2010], w=all_atus$weight_adj[all_atus$year==2010])
m11 <- mean(all_atus$home_paper[all_atus$year==2011], w=all_atus$weight_adj[all_atus$year==2011])
m12 <- mean(all_atus$home_paper[all_atus$year==2012], w=all_atus$weight_adj[all_atus$year==2012])
m13 <- mean(all_atus$home_paper[all_atus$year==2013], w=all_atus$weight_adj[all_atus$year==2013])
m14 <- mean(all_atus$home_paper[all_atus$year==2014], w=all_atus$weight_adj[all_atus$year==2014])
mean <- cbind(m03,m04, m05, m06, m07, m08, m09, m10, m11,m12, m13, m14)

year <- seq(2003, 2014, 1)

#Format properly
mat2<- rbind(year, mean)
mat2 <- t(mat2)
colnames(mat2)[colnames(mat2)==""] <- "mean"
mat2.d <- as.data.frame(mat2)

#plot for nonmarket work ####
ggplot(mat2.d, aes(x=mat2[,1], y=mat2[,2])) +
  geom_point(shape=1) + 
  geom_smooth(method=lm, linetype = 2, se= FALSE) +
  geom_line(size = 1) +
  #       scale_y_continuous(limits = c(20, 30)) +
  scale_x_continuous(breaks = 2003:2014) +
  xlab("Year") +
  ylab("Hours Per Week") +
  ggtitle("Nonmarket Work")



#For lesiure: 
#Take the mean of every year, and create a matrix for plot:
m03 <- mean(all_atus$leisure_paper[all_atus$year==2003], w=all_atus$weight_adj[all_atus$year==2003])
m04 <- mean(all_atus$leisure_paper[all_atus$year==2004], w=all_atus$weight_adj[all_atus$year==2004])
m05 <- mean(all_atus$leisure_paper[all_atus$year==2005], w=all_atus$weight_adj[all_atus$year==2005])
m06 <- mean(all_atus$leisure_paper[all_atus$year==2006], w=all_atus$weight_adj[all_atus$year==2006])
m07 <- mean(all_atus$leisure_paper[all_atus$year==2007], w=all_atus$weight_adj[all_atus$year==2007])
m08 <- mean(all_atus$leisure_paper[all_atus$year==2008], w=all_atus$weight_adj[all_atus$year==2008])
m09 <- mean(all_atus$leisure_paper[all_atus$year==2009], w=all_atus$weight_adj[all_atus$year==2009])
m10 <- mean(all_atus$leisure_paper[all_atus$year==2010], w=all_atus$weight_adj[all_atus$year==2010])
m11 <- mean(all_atus$leisure_paper[all_atus$year==2011], w=all_atus$weight_adj[all_atus$year==2011])
m12 <- mean(all_atus$leisure_paper[all_atus$year==2012], w=all_atus$weight_adj[all_atus$year==2012])
m13 <- mean(all_atus$leisure_paper[all_atus$year==2013], w=all_atus$weight_adj[all_atus$year==2013])
m14 <- mean(all_atus$leisure_paper[all_atus$year==2014], w=all_atus$weight_adj[all_atus$year==2014])
mean <- cbind(m03,m04, m05, m06, m07, m08, m09, m10, m11,m12, m13, m14)

year <- seq(2003, 2014, 1)

#Format properly
mat3<- rbind(year, mean)
mat3 <- t(mat3)
colnames(mat3)[colnames(mat3)==""] <- "mean"

mat3.d <- as.data.frame(mat3)

#plot for leisure
#plot for nonmarket work ####
ggplot(mat3.d, aes(x=mat3[,1], y=mat3[,2])) +
  geom_point(shape=1) + 
  geom_smooth(method=lm, linetype = 2, se= FALSE) +
  geom_line(size = 1) +
  #       scale_y_continuous(limits = c(20, 30)) +
  scale_x_continuous(breaks = 2003:2014) +
  xlab("Year") +
  ylab("Hours Per Week") +
  ggtitle("Leisure")


#Going forwards: use mat in plot and then put trend over top that.e.g.:
plot(mat[,1], mat[,2], type = "o")


########Column 6s############
#Run Table 1, column 3 right before this.
#Table 1, column 6: 
new <- data.frame(year = c(2013, 2014))

a<- predict(trend1, newdata= new)  #Get predicted points
a <-sum(a)/2 #average 2013/14
a<- a-out1$coefficients ###This is the result
a


trendb<- lm(worka_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
b<- predict(trendb, newdata= new) 
b <-sum(b)/2 
b<- b-out2$coefficients 
b

trendc<- lm(worku_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
c<- predict(trendc, newdata= new) 
c <-sum(c)/2 
c<- c-out3$coefficients 
c

trendd<- lm(childcare_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
d<- predict(trendd, newdata= new) 
d <-sum(d)/2 
d<- d-out4$coefficients 
d

trende<- lm(home_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
e<- predict(trende, newdata= new) 
e <-sum(e)/2 
e<- e-out5$coefficients 
e

trendf<- lm(homeproduction_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
f<- predict(trendf, newdata= new) 
f <-sum(f)/2 
f<- f-out6$coefficients 
f

trendg<- lm(homeown_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
g<- predict(trendg, newdata= new) 
g <-sum(g)/2 
g<- g-out7$coefficients 

trendh<- lm(shopping_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
h<- predict(trendh, newdata= new) 
h <-sum(h)/2 
h<- h-out8$coefficients 

trendi<- lm(othercare_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
i<- predict(trendi, newdata= new) 
i <-sum(i)/2 
i<- i-out9$coefficients 

trendj<- lm(leisure_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
j<- predict(trendj, newdata= new) 
j <-sum(j)/2 
j<- j-out10$coefficients 

trendk<- lm(tv_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
k<- predict(trendk, newdata= new) 
k <-sum(k)/2 
k<- k-out11$coefficients

trendl<- lm(socializing_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
l<- predict(trendl, newdata= new) 
l <-sum(l)/2 
l<- l-out12$coefficients 

trendm<- lm(sleeping_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
m<- predict(trendm, newdata= new) 
m <-sum(m)/2 
m<- m-out13$coefficients 

trendn<- lm(ep_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
n<- predict(trendn, newdata= new) 
n <-sum(n)/2 
n<- n-out14$coefficients 

trendo<- lm(otherleisure_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
o<- predict(trendo, newdata= new) 
o <-sum(o)/2 
o<- o-out15$coefficients 

trendp<- lm(other_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
p<- predict(trendp, newdata= new) 
p <-sum(p)/2 
p<- p-out16$coefficients 

trendq<- lm(education_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
q<- predict(trendq, newdata= new) 
q <-sum(q)/2 
q<- q-out17$coefficients 

trendr<- lm(civic_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
r<- predict(trendr, newdata= new) 
r <-sum(r)/2 
r<- r-out18$coefficients 

trends<- lm(ownmedical_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
s<- predict(trends, newdata= new) 
s <-sum(s)/2 
s<- s-out19$coefficients 


#Table 2, column 6:
#Do not run code aove after code below is run.
#need to rerun code to get new data
#and trends1 through 3
all_atus <- all_atus[all_atus$male==1,]
####NOTE: IT IS VERY IMPORTANT to run Table 2, column 3 before running the code below
#alright

trenda<- lm(work_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
a<- predict(trend1, newdata= new)  #Get predicted points
a <-sum(a)/2 #average 2013/14
a<- a-out1$coefficients ###This is the result
a

trendb<- lm(worka_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
b<- predict(trendb, newdata= new) 
b <-sum(b)/2 
b<- b-out2$coefficients 

trendc<- lm(worku_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
c<- predict(trendc, newdata= new) 
c <-sum(c)/2 
c<- c-out3$coefficients 

trendd<- lm(childcare_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
d<- predict(trendd, newdata= new) 
d <-sum(d)/2 
d<- d-out4$coefficients 

trende<- lm(home_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
e<- predict(trende, newdata= new) 
e <-sum(e)/2 
e<- e-out5$coefficients 

trendf<- lm(homeproduction_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
f<- predict(trendf, newdata= new) 
f <-sum(f)/2 
f<- f-out6$coefficients 

trendg<- lm(homeown_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
g<- predict(trendg, newdata= new) 
g <-sum(g)/2 
g<- g-out7$coefficients 

trendh<- lm(shopping_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
h<- predict(trendh, newdata= new) 
h <-sum(h)/2 
h<- h-out8$coefficients 

trendi<- lm(othercare_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
i<- predict(trendi, newdata= new) 
i <-sum(i)/2 
i<- i-out9$coefficients 

trendj<- lm(leisure_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
j<- predict(trendj, newdata= new) 
j <-sum(j)/2 
j<- j-out10$coefficients 

trendk<- lm(tv_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
k<- predict(trendk, newdata= new) 
k <-sum(k)/2 
k<- k-out11$coefficients

trendl<- lm(socializing_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
l<- predict(trendl, newdata= new) 
l <-sum(l)/2 
l<- l-out12$coefficients 

trendm<- lm(sleeping_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
m<- predict(trendm, newdata= new) 
m <-sum(m)/2 
m<- m-out13$coefficients 

trendn<- lm(ep_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
n<- predict(trendn, newdata= new) 
n <-sum(n)/2 
n<- n-out14$coefficients 

trendo<- lm(otherleisure_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
o<- predict(trendo, newdata= new) 
o <-sum(o)/2 
o<- o-out15$coefficients 

trendp<- lm(other_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
p<- predict(trendp, newdata= new) 
p <-sum(p)/2 
p<- p-out16$coefficients 

trendq<- lm(education_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
q<- predict(trendq, newdata= new) 
q <-sum(q)/2 
q<- q-out17$coefficients 

trendr<- lm(civic_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
r<- predict(trendr, newdata= new) 
r <-sum(r)/2 
r<- r-out18$coefficients 

trends<- lm(ownmedical_paper ~ year, data=subset(all_atus, year<=2008), weights=weight_adj) 
s<- predict(trends, newdata= new) 
s <-sum(s)/2 
s<- s-out19$coefficients 

